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Abstract 

We investigate a nonlinear coherent feedback circuit constructed from pre-existing supercon- 
ducting microwave devices. The network exhibits emergent bistable and astable states, and we 
demonstrate its operation as a latch and the frequency locking of its oscillations. While the net- 
work is tedious to model by hand, our observations agree quite well with the semiclassical dynamical 



model produced by a new software package [N. Tezak et ai, arXiv:1111.3081vl] that systematically 



interpreted an idealized schematic of the system as a quantum optic feedback network. 
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The degree of control over matter and electromagnetic fields demonstrated in the past two 
decades suggests that quantum engineering may become a powerful discipline. However, the 
extreme requirements for quantum scale engineering, be it for quantum- pQ or ultra-low 
energy classical-information systems [2], suggest that active feedback will be necessary in 
useful networks [H EJ Hj. But while important proof-of-principle demonstrations of quan- 
tum error correction (QEC) have been reported for instance [5H7], the unwieldy classical 
feedback equipment so far employed poses perhaps the greatest obstacle to realizing useful, 
complex systems. More generally, measurement-based quantum feedback [8HT2"] may prove 
impractical simply because "measurement" implies a network interruption by a fundamen- 
tally non-integrable system. To overcome this bottleneck, quantum networks may need 
to actively stabilize themselves through coherent feedback of probes without measurement 
El [T3H20] . Moreover, [HI [15] suggests that coherent feedback can outperform even ideal 
measurement-based feedback. 

Measurement-based feedback to superconducting microwave quantum circuits is partic- 
ularly difficult as signal transfer between a cryostat and room temperature electronics is 
inefficient and slow (TTJ [121 EI]- It was proposed in [IB] that coherent feedback circuits 
employing non-linear (Kerr) resonators are a natural approach to self-stabilizing, digital 
optical information processing in a complex quantum network. Here, we demonstrate that 
these insights readily apply to superconducting circuits by constructing a coherent feed- 
back multivibrator network (a circuit operable as a set-reset latch or an astable oscillator) 
from pre-existing Kerr-type resonators and coherent feedback of signals that never leave the 
<50mK environment. This network becomes useful when integrated with other systems, 
and could act as a binary controller in a larger QEC coherent feedback network [IHl H7] or 
as a cryogenic clock. And while an idealized model of this device could be derived manually, 
more complex systems would prove intractable. Thus, we demonstrate that our observations 
agree quite well with a semiclassical model that was systematically produced from a network 
schematic by a hierarchical quantum circuit modeling package [22]. While previous experi- 
ments have validated similar approaches to modeling coherent feedback circuits in linear [19] 
and linear-quantum [20] optical networks, to our knowledge this is the first application to a 
nonlinear network, in a superconducting microwave context, and using automated quantum 
circuit modeling. 

The network's primary components are two single port microwave resonant circuits whose 
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resonance frequency is power dependent and tunable with an applied magnetic flux 
These tunable Kerr circuits (TKCs) were originally fabricated to serve as Josephson para- 
metric amplifiers for near quantum-limited amplification of weak microwave signals and the 
preparation of squeezed microwave fields [25], [26] (see also [27]). The TKCs are quarter-wave 
transmission line resonators formed by a coplanar waveguide with one end shorted and a 
capacitively coupled port at the other, and were mounted in separate sample boxes. A series 
array of 40 Josephson junction SQUIDS interrupt the coplanar waveguide center conductor, 
providing a non-linearity that makes the devices' input-output (I/O) properties analogous 
to that of a high-quality, single-sided optical Kerr cavity (with Kerr coefficient x < 0) [23J . 
Thus, the reflected phase is a non-linear function of input power [24] . This function can 
even be bistable for input drives that simultaneously are detuned below the TKCs' center 
frequency coo by at least the critical value ojq — co PtC = A c = a/3k and exceed the critical power 
P c = fthjp x 4/t 2 /(3v / 3|xD? where k is the field decay rate of the TKC [21]. The TKCs used 
here both have k/2h = 15 MHz and P c = —98 ± 2 dBm (uncertainty in the line calibration) 
when tuned such that u /2ii = 6.408 GHz. 

When the input drive detuning is close to, but does not exceed A c , a TKC is monostable 
for all input powers, but the phase of the reflected signal 'flops' by approximately ir radians 
when the power, p x P c , exceeds P c , see Fig. [pi. Because such phase shifts are readily 
converted into power variations in an interferometric network, [18] suggested that Kerr 
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FIG. 1: (color online) a) Each TKC operates as a single-sided Kerr cavity. Reflected phase is a 
nonlinear function of drive amplitude for drive detunings < A c . b) Network schematic. Two flux 
biased TKCs are connected with approximately 25 cm cable connections to a quadrature hybrid 
in a feedback configuration. Phase locked signals drive the ports Ino and Ini, and are separated 
from output signals by directional couplers. 



3 



cavities may usefully approximate NAND gates in an optical network. For example, if two 
channels carrying either p > 1/4 (high) or p < 1/4 (low) interfere in phase and are directed 
at a Kerr cavity, the phase of the reflected signal will 'flop' only if both inputs are high. 
Moreover, while a NAND gate/Kerr cavity is monostable in isolation, a network of two 
NAND gates/Kerr cavities in mutual feedback may function as a multivibrator. Adapting 
[18] to our context, a coherent feedback network of two TKCs should display emergent 
bistable and astable dynamics. Such a network would also be nearly lossless and suitable 
for chip-level integration with quantum information microwave systems. 

Represented in Fig. [TJd, the network components are housed in a dilution refrigerator and 
consist of two TKCs and a 4-8 GHz commercial quadrature hybrid (analogous to an optical 
50/50 beamsplitter). The TKCs are connected to the hybrid in such a way that signals 
they reflect are split between one of the two network outputs and the other TKCs input, 
producing a coherent feedback network. These connections were made by low-loss, coaxial 
Cu cables, but our lab has previously interconnected these components on a single chip 
[25] . Two signal generators drive the system through low-temperature attenuation stages, 
producing two phase—locked, low-temperature microwave drive inputs to the network. The 
signals reflected out these same lines are separated from the inputs by directional couplers, 
and are amplified by two low-noise cryogenic HEMTs for analysis. 

Superconducting microwave devices are often describable with models equivalent to I/O 
models in quantum optics [29, 30J. In such cases (e.g. TKCs and hybrids), one may model 
interconnected devices using cascaded I/O techniques still developing in quantum optics JT3J, 
I3T1 132] . Unfortunately, these calculations are tedious, even for networks as basic as Fig. [T|d. 
A new software package, Quantum Hardware Description Language (QHDL) [22], adapts a 
standard electrical engineering modeling language to automate this modeling, interpreting a 
schematic diagram input that specifies the bosonic field I/O connections between pre-defined 
quantum optical primitive or composite models. 

Here, after a schematic representing Fig. [TJd is loaded, QHDL outputs the network's 
symbolic Heisenberg equations of motion (EOM). As these TKCs typically contain ~1000 
photons, driven by coherent fields, a semiclassical approximation is invoked to simulate mean 
field dynamics, normal ordering the system operators in the EOM and replacing operators 
with complex scalars [33] ■ While these approximations could have been applied at the device 
level, the quantum network model construction is no more difficult than the coherent classical 
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one. In our approach, future work considering the effects of intrinsic quantum fluctuations 
or integration with necessarily quantum systems follows readily. QHDL employs a standard 
approach to I/O theory that assumes that transmission line delays are negligible [T3 t I5T | 152] . 
Furthermore, to compare our experiment to the ideal, in our model the TKCs are lossless 
and identical single-sided Kerr cavities, interconnections produce symmetric phase shifts 
and loss, and unwanted reflections are ignored. 

In testing our model, we first measure the network's linear behavior by probing it with 
power p < 1. Tuning the TKCs far outside the probed region, we observe interferometric 
resonances, with the input power periodically distributed between the two outputs as the 
drive frequency varies. With the TKCs tuned to be co-resonant near the middle of the 
probe range, additional resonances appear, and avoided crossings between all resonances 
are apparent, indicating that the TKCs are coherently coupled at rate ~ k to each other 
and to the network (Fig[2^-b). Comparing the data with model simulations (Fig. j^-d), 
we calibrate 0.4 dB round-trip loss in each interconnection. We note that the network's 
interconnections are longer than needed, a compromise between wanting long connections 
so that a desired phase shift between components could be achieved through frequency 
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FIG. 2: (color online) Magnitude and phase response data (a & b) and simulation (c & d) in 
the linear, p <C 1 regime. Driving Ini only and with the TKCs out of probe range, the response 
measured at Outn (Outi) is shown in purple (red). Gold (cyan) lines depict Outo (Outi) responses 
when both TKCs are co-resonant at oj/2-k = 6.408 GHz (dotted vertical lines). 
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tuning and wanting short connections such that the delay between components be negligible. 
Intending to consider dynamics only on time scales greater than when the experiment 
was deployed, we chose 25 cm interconnections (resulting in a .24k -1 delay between TKCs), 
producing a 385 MHz period in the frequency response. 

Beyond the linear regime, and fixing both TKCs to ujq/2tt = 6.408 GHz and our probe 
frequency to u p /2n = 6.39 GHz (so that uq — uj p = 0.69A C and the round trip phase shift per 
cable is 2.65 rad), the network exhibits both bistable and astable dynamics. For balanced, 
p > 1 and out-of-phase drives on the two inputs, coherent feedback causes the network to 
be bistable, with either high power driving TKCo and Outi and low power driving TKCi 
and Outo, or vice versa. To give a heuristic explanation (see [33] for a quantitive model), 
if the power incident on TKCo happens to be high (p > 1), it is reflected with a 'flopped' 
phase shift (~ tt rad). For these biasing conditions, the TKCo-reflected signal interferes with 
the Inx drive at the hybrid such that more power is directed to Outi than TKCi. The low 
power (p < 1) signal incident on TKCi is then reflected with no additional phase shift, and 




FIG. 3: (color online) a-b) Mean Outo power, adiabatically sweeping either Ino's or Ini's drive 
amplitude high-and-low, with different biasings on the other input, c) Alternately sweeping in 
both directions produces a colored mesh (color indicating Outo power) where the bistable region 
is apparent where different colored lines intersect. White lines demarcate the observed bistability 
region; brown dot represents the 'hold' biasing state; arrows represent 'set' or 'reset' operations, 
d) Simulation of c. 
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consequently interferes at the hybrid with the Ino drive such that more power drives TKCo 
than Outo, reinforcing the original, strong TKCo drive. By symmetry, for the same biasing 
conditions, the opposite network state is also self-stabilizing. Thus, while both TKCs would 
be monostable in isolation at this detuning, the network exhibits a bistable output regime 
when the two input drives are balanced and strong. If one further increases the Ino (I n i) 
drive enough relative to the other, bistability disappears, and the system relaxes to a high 
Outx (Out ) state. 

In Fig. |3^L-b, we plot the mean Out power observed as a function of the input drives, 
as the amplitude of either the In or In! drive is adiabatically swept high-and-low at 1 kHz 
for 100 cycles while the other input is fixed at various amplitudes. Several hysteresis loops 
are apparent in the regime where the the two input drives are roughly equal, a consequence 
of the bistable dynamics described above. The Outi powers are largely symmetric upon 
exchange of the input axes (asymmetries being a consequence of slight network asymmetries 
not considered here). The bistable region as a function of the two inputs becomes clearer 
in Fig. ^jp, where low-to-high and high-to-low sweeps of one input are alternated for various 
static biasings of the other input, and the mean Outo power is depicted on the same color 
scale. This produces a colored mesh that indicates the bistable region by the intersection 
of different color lines. Fig. [3jd is a simulation of the same, producing a very similar color 
pattern and a similar, but more symmetric bistability region. All output power data was 
calibrated first by scaling the signal measured at Outi such that for far-detuned TKCs and 
balanced inputs (inferred by the 6.3-6.55 GHz phase response) the output powers were equal 
(compensating for amplifier asymmetries), then by equally scaling both outputs such that 
the highest Outo power in Fig. matched the highest simulated power in Fig. |3ji. 

This bistability may be leveraged to operate the network as a set-reset latch (or 'flip- 
flop'), a binary memory element that outputs power according to prior inputs [18]. In Fig. [4^l, 
the averaged output response is tracked as the two input drives are amplitude modulated 
between a 'hold' condition of equal, p = 1.6 drives, and either the Ino ('set') or Ini ('reset') 
drives doubling in power and returning. Fig. |3Jd simulates the same. The hold condition 
corresponds to the brown dot in Fig. ^jp, while the set and reset operations correspond to 
modulating the input powers according to the horizontal and vertical arrows, respectively. 
As the hold state is bistable and connected to the monostable set and reset states via different 
stable manifolds (Fig. [3]), each set-hold (reset-hold) event causes the Outi (Outo) signal to 
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swing high regardless of the prior state. While the modulation frequency is again ~kHz, 
the network's response rate is at least that of the 2 MHz detection bandwidth. To note 
one potential application, [TBI HZ] suggests that set-reset sub-networks like these could act 
as binary controllers in 'hard-wired' implementations of QEC, stabilizing superconducting 
qubit arrays in a larger coherent feedback network. 

Increasing the detection bandwidth to 50 MHz, various drive settings produce sustained 
output power oscillations at frequencies ~ k. For example, Fig. [5^, represents the mixed- 
down power spectrum detected at Out while driving only In with a continuous wave 6.39 
GHz tone of various amplitudes. Starting near p = —1 dB, ~10 MHz and higher harmonics 
emerge and accelerate with input power. To compare to the bistable case, with only one 
drive, a strong or weak signal reflected by TKCo has no drive to interfere with. Thus, as 
TKCo equilibrates, TKCi is driven with a relatively strong or weak signal, respectively, 
the opposite of the bistable case. Consequently, when this signal ultimately reflects back 
towards TKCo, the network destabilizes and oscillates between both states [33J. 

In this case, however, the analogous simulation (Fig. predicts 17 MHz power oscilla- 
tions, their emergence at p — 1.2 dB, and their increasing frequency with drive power. In 
view of the accuracy of the idealized simulations to reproduce the low-frequency dynamics 
in Figs. [l]|4]and the significant frequency-dependence of phase shifts of dynamic signals >10 
MHz in our physically-extended network (see Fig. |2|, we suspect the discrepancy stems from 
the zero delay assumption of QHDL. This hypothesis is supported by a linearized version of 
the QHDL-derived model. Analysis of the dynamics about the EOM fixed points predicts 
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FIG. 4: (color online) Mean output response data (a) and simulation (b) depicting the network's 
operation as a binary memory element. After each 'set'/'reset' operation (blue S/red R), input 
power is primarily directed out Outi/Outn, even after retuning to the 'hold' state (green H). Input 
states of variable durations were an experimental convenience. 
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the emergence of a stable 17 MHz limit cycle at 1.2 dB, exactly as observed in the Fig. 
simulation. Adding the approximate effects of transmission line delays to the linearized 
model destabilizes the dynamics at lower drive powers, suggesting a stable 11 MHz limit 
cycle emerging at -1.6 dB [33], much closer to what is experimentally observed (Fig. [5^,). I/O 
models may be generalized to include finite delays and while the resulting models may be 
automated, they were deemed too complex for first generation software [131 E2] . Cascaded 
I/O models are most appropriate for chip-scale systems as opposed to our extended net- 
work; chip-scale integration would improve simulation accuracy if our hypothesis is correct. 
It is worth mentioning, though, that QHDL's qualitative accuracy beyond its range of strict 
applicability was quite useful for predicting astable parameter regions. 

Finally, we demonstrate (measurement-based) stabilization of these oscillations in 
Fig. [5ji. By setting the Ino drive to p — 2.58 dB and mixing down the Outi signal with a 6.4 
GHz local oscillator (10 MHz detuned from the injected tone) significant phase noise relative 
to our room temperature frequency standard is apparent (likely due to technical jitter in 



Out. power spectrum data c ) Out power spectrum simulation 

50 — 




10 

Frequency [MHz] 



FIG. 5: (color online) a) Mixed-down power spectrum detected at Outn as Ino is driven at 6.39 
GHz, input powers relative to P c . b) Outn power oscillations in time with 2.58 dB input, c) 
Simulation of the Outn power spectrum. Red (blue) circle marks the predicted emergence of a 
stable limit cycle when delays are not (are) added to a linearized model [33J. d) Power spectrum 
of the Outn signal frequency component 10 MHz detuned from the injected tone with and without 
frequency locking. 



9 



the TKC center frequencies). As the frequency of the output power oscillations varies with 
input power, using this phase signal to drive .23 dB analog amplitude modulation of the 
injected tone (100 kHz modulation bandwidth) creates a phase locked loop that stabilizes 
the 10 MHz pulse train spontaneously produced by our cryogenic network to the 10 MHz 
room temperature clock that phase locks our generators. 

While these dynamics are classical, QHDL outputs quantum models and TKCs are rou- 
tinely used by our lab to generate and measure non-trivial quantum fields [25]. It would 
be interesting, for instance, to consider how quantum field fluctuations propagate through 
this network and perhaps disturb the mean-field dynamics reported here [34] . Nonethe- 
less, classical dynamics are sufficient to demonstrate that classical information systems are 
readily produced by coherent feedback on generic quantum devices. But because they are 
constructed from the same hardware as quantum microwave circuits, they hold a natural 
advantage in terms of the chip-level classical/quantum integration that would be necessary 
for truly scalable quantum circuits [TBI E] • We conclude by reiterating that this system was 
constructed from pre-existing components of types generically available in superconducting 
circuit labs [231 EZ] • And while this system's intricate and potentially useful dynamics are 
difficult to consider manually, they are readily analyzed and integrated into larger network 
models using a laptop and a small number of I/O laws originally formulated for quantum 
optics. This observation suggests that automated modeling techniques like QHDL are now 
needed to properly compliment quantum hardware advances. 

We acknowledge partial support from the DARPA QuEST program and from the NSF 
Physics Frontier Center. JK acknowledges the NRC for financial support, W. Kindle and 
H.-S. Ku for experimental advice, and N. Tezak and H. Mabuchi for very helpful discussions 
and the beta version of QHDL. 



* Electronic address: jkerc@jila.colorado.edu 
[1] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quanutm Information (Cam- 
bridge University Press, 2000). 
[2] D. A. B. Miller, Nat. Photon. 4, 3-5 (2010). 
[3] N. C. Jones et al, ||arXiv:1010.5022[ 



10 



[4] H. Mabuchi, Appl. Phys. Lett. 98, 193109 (2011). 

[5] M. D. Reed et al, Nature 482, 382-385 (2012). 

[6] P. Schindler et al, Science 332, 1059 (2011). 

[7] T. Aoki et al, Nat. Physics. 5, 541-546 (2009). 

[8] W. P. Smith et al, Phys. Rev. Lett. 89, 133601 (2002). 

[9] M. A. Armen et al, Phys. Rev. Lett. 89, 133602 (2002). 

[10] C. Sayrin et al, Nature 477, 73-77 (2011). 

[11] R. Vijay et al, |arXiv:1205.5591[ 



[12] D. Riste, C. C. Bultink, K. W. Lehnert, and L. DiCarlo, [arXiv: 1207.2944] 



[13] J. Gough and M. R. James, IEEE Trans. Auto. Control. 54, 2530-2544 (2009); Comm. in 

Math. Phys. 287, 1109-1132 (2008). 
[14] H. I. Nurdin, M. R. James, and I. R. Petersen, Automatica 45, 1847 (2009). 



[15] R. Hamerly and H. Mabuchi, arXiv:1206.2688 arXiv:1206.0829 



[16] J. Kerckhoff, H. I. Nurdin, D. S. Pavlichin, and H. Mabuchi, Phys. Rev. Lett. 105, 040502 
(2010). 

[17] J. Kerckhoff, D. S. Pavlichin, H. Chalabi, and H. Mabuchi, New J. Phys. 13 055022 (2011). 
[18] H. Mabuchi, Appl. Phys. Lett. 99, 153103 (2011). 
[19] H. Mabuchi, Phys. Rev. A 78, 032323 (2008). 

[20] S. Iida et al, IEEE Trans. Auto. Control, (accepted) (2012); [arXiv: 1103. 1324^ 1. 
[21] I. Siddiqi, Supercond. Sci. Technol. 24, 091002 (2011). 



[22] N. Tezak et al, Phil. Trans. Roy. Soc. A (accepted) (2012); |arXiv: 1111 .3081^ 1. 
[23] M. A. Castellanos-Beltran et al, Nat. Phys. 4, 929 (2008). 
[24] B. Yurke and E. Buks, J. Lightwave Technol. 24, 5054 (2006). 
[25] F. Mallet et al, Phys. Rev. Lett 106, 220502 (2011). 
[26] J. D. Teufel et al, Nature 475, 359363 (2011). 

[27] B. Abdo et al, Appl. Phys. Lett. 99, 162506 (2011); R. Vijay, M. H. Devoret, and I. Siddiqi, 
Rev. Sci Instrum. 80, 111101 (2009); R. Vijay, D. H. Slichter and I. Siddiqi, Phys. Rev. 
Lett. 106, 110502 (2011). 

[28] H. S. Ku et al, IEEE Trans. Appl. Superconductivity 21, 452-455 (2011). 

[29] B. Yurke and J. S. Denker, Phys. Rev. A 29, 1419 (1984). 

[30] A. A. Clerk et al, Rev. Mod. Phys. 82, 1155-1208 (2010). 

11 



[31] C. W. Gardiner and P. Zoller, Quantum Noise (Springer, 2004). 

[32] H. J. Carmichael, Phys. Rev. Lett. 70, 2273-2276 (1993). 

[33] See Supplemental Material at [URL will be inserted by publisher] for modeling details. 

[34] J. Kerckhoff, M. A. Armen, and H. Mabuchi, Opt. Express 19, 24468-24482 (2011). 

[35] C. W. Gardiner, Phys. Rev. Lett. 70, 2269-2272 (1993). 

[36] S. H. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chem- 
istry and Engineering (Perseus, Cambridge, MA, 1994). 

[37] C. Jeffries and K. Wiesenfeld, Phys. Rev. A 31, 1077 (1985); Phys. Rev. A 33, 629 (1986). 

SUPPLEMENTARY INFORMATION 
General modeling 

The power of the Quantum Hardware Description Language (QHDL) [22] modeling ap- 
proach (which automates the quantum circuit algebra of Gough and James [13], which in 
turn generalizes earlier work on cascaded open quantum systems by Carmichael [32] and 
Gardiner [3"T| 133] ) stems from the fact that individual open quantum optical components 
are given the same succinct representation as interconnected networks of quantum optical 
components. 

This representation consists of a triple (S,L,if). To describe briefly, S is an operator- 
valued, square scattering matrix that specifies how input (uni-directional), freely-propagating 
bosonic fields are directly scattered to output fields (as many vector indices as there are 
input fields), as in the action of a beamsplitter or microwave hybrid. L is an operator- 
valued coupling vector that specifies how each field mode couples to the internal quantum 
degrees of freedom (if any) of component devices. And H is the effective Hamiltonian that 
specifies the internal dynamics of the devices, independent of the effects of the free fields. 

The construction of a network (S, L, H) from component triples proceeds with a small set 
of composition rules (here presented assuming negligible time delay between components), 
depicted in Fig. [6] The concatenation product represents the effective dynamics of two 
components that have no direct free field interconnection, but could share a common internal 
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Hilbert space: 

(SlH2, LlBB2> Hia2) = (Si,Li,i?i) ffl (S2, L2, #2) 





Si 











( 













s 2 




^2 





(1) 



The series product represents the effective dynamics of a network in which the output fields 
of component 2 are fed into the inputs of component 1 



vSi<, 2 , L 1<2 , H\<Q t 



(S 1? L 1; H x ) < (S 2) L 2 , F 2 ) = (SxSs, L x + S^, H 2 + H x + ^LjSiLa} 



(2) 

where ^{^4} = (A — A*)/2i and the ' operation returns a transposed operator matrix 
with operator adjoints in its entries. Finally, the feedback operation represents the effective 
network dynamics when the k th output channel is fed back into the I th input (thus reducing 





N m2 = N x ffl N 2 N 1<2 = N 1 <N 2 
a) b) 







FIG. 6: Depictions of the essential composition operations through which component representa- 
tions are combined to form composite network representations, a) The concatenation product, b) 
The series product, c) The feedback operation. Adapted from 
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the number of input and output ports by 1): [(S, L, H)] 



k-tl 



(S, L, H, ) where 



Si,i 

Sfe-1,2 

Sk+i,i 

S n ,i 
Si,i 

Sk-1,1 
Sk+i,i 



:i-s, 



kdj 



Sk,i ■ ■ ■ Sk,i-i Sk,i+i ■ ■ ■ Sk,i 



s. 



n,l 



>1 — Ski) 1 Lk 



H = H + % 



.3=1 



(1 — Sk,i) x Lk 



(3) 



tii 



where Sj^g- and indicate the original scattering matrix and coupling vector with the k 
row and I th column removed. For more details of the fundamental models and assumptions, 
we refer readers to [13, 22J. 

Whether a (S, L, H) triple describes an individual component or a network of components, 
the effective dynamics of the system are calculated in the same way For example, assuming 
the input fields are in the vacuum state, the evolution of an operator X that acts on the 
internal Hilbert space (e.g. the annihilation operator of a TKC mode) is systematically 
calculated as [13] (h = 1) 



dX 



l[X,H] + ^V[X,L] + ^[L\X]L^j dt+dA^(t)S^[X,L)+[V,X]SdA(t)+Ti [(S f XS - X)dA T (t)] 

(4) 



where T is the operator matrix transpose. A(t), A^(t) are operator vectors, whose entries are 
known as as quantum noise processes, whose infinitesimal increments (e.g. dA[k](t)) may be 
roughly considered the annihilation and creation operators (respectively) on the infinitesimal 
segment of input free field that interacts with the component or network at time t. A is an 
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operator matrix whose entries are a third kind of quantum noise process whose increments 
may be roughly considered bilinear products of field annihilation and creation operators (e.g. 
its diagonal elements are similar to number operators on each infinitesimal field segment). 
Also, the output fields are related to the input fields and the internal degrees of freedom by 

dA out (t) = SdA(t) + Ldt, (5) 

as well as related relations for A*(t) and A(t). 

Thus, when the assumptions are valid, the dynamics of both individual quantum optical 
components and complex networks of interconnected components may be derived systemati- 
cally: following a schematic of interconnected (S, L, H) models, one first derives the effective 
(S, L, H) for the entire network using rules Eqs. (??-[3]); then, one derives the quantum equa- 
tions of motion using Eqs. (??-??). Often, however, this general procedure is very tedious. 

The most immediate value of the Quantum Hardware Description Language (QHDL) 
[22] is that it insulates a user from this computational tedium. One may produce the desired 
equations of motion from an intuitive schematic diagram and less than 10 lines of code. 

Specific model 

Following this general modeling and using procedures analogous to [Ml ESI ED], one may 
derive the T = (Stkc, ^tkc, Htkc) triple representation for an ideal TKC as a single mode 
component: 

$tkc = [—1] 
Ltkc = [— iV^Ka] 

H TKC = Aa f a + |a f V (6) 

where a is the annihilation operator on the TKC resonator mode, A = ojq — u p is the 
detuning between the TKC resonance frequency (ojq) and the carrier frequency of the input 
field driving the TKC (oo p ), k is the field decay rate, and x < is the effective Kerr coefficient 
produced by the SQUID array. The remaining component types employed in the network 
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model are: beamsplitters BS = (Sbs,^bs, Hbs) 



v fi 



J BS 



where |/i| 2 + |z/| 2 = 1; phase shifters $ = (S^,, L^, if^) 



H 



BS 



(7) 



S = [e^], L =[O], ff* = 0; 



(8) 



and coherent drives W a = (Swa> Wai H\ 



Wai 



&Wa — [1], LvKa — [«]> #Wcn — 0, 



(9) 



with complex amplitude a. From these general component models, the TKCs are taken as 
distinct but identical T components. The quadrature hybrid is modeled as the concatenation 
of two beamsplitters, H ~ BSqSBSi, with appropriate relations between the reflection and 
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FIG. 7: Schematic representation of the network model employed in the main Letter and in- 
terpreted by QHDL. Overall a 4-mode input-output network (2 signal channels, 2 loss channels), 
individual components are icons representing quantum optical (S,L,i7) models, with connections 
between components representing (uni-directional) bosonic field modes. Coherent drive 'compo- 
nents' are not shown, but are eventually placed upstream of Ino and Ini ports. This schematic and 
the sub-componets it references define the network model returned by QHDL. 
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transmission coefficients {/iq, uq} and {[ii,Vx} stemming from the fact that this single, bi- 
directional physical component is modeled as two uni-directional beamsplitters (with the 
"~" representing the fact that some field index re-ordering is also employed). Transmission 
line-induced phase shifts are modeled as two identical $ components, and transmission line 
loss is modeled by two identical beamsplitters that mix the transmission line modes with 
vacuum at a low rate, i.e. \u\ <C The two coherent drives are modeled as two W a 
components 'upstream' of the network, which displace input vacuum fields by respective 
amplitudes. 

Icons that represent these components are arranged in the schematic diagram shown in 
Fig. [7j with interconnections that emulate our experimental network. From this schematic, 
QHDL parsers were employed to calculate first the effective (S, L, H) representation of the 
network and then semiclassical approximations of its equations of motion. To give a concrete 
example of the calculation procedure, we will devote most of the remainder of this section to 
outlining the procedure and results obtained in the case of an even simpler network model 
that is lossless and has integer 7r-radian phase shifts. 

If one removes the 'LossO' and 'Lossl' beamsplitter components and associated input and 
output ports in the Fig. [7] schematic, the network model without any coherent drives may 
be characterized as 



where we have introduced two new types of (S, L, H) 'components' necessary for appropriate 
field indexing: the permutation matrix -P(i,o) that reverses the ordering of the two output 
fields and the identity component I n that passes n-input modes to outputs without scaling 
or re-ordering. In plain English this sequence may be read as 

"Output 4 of if is fed into $1 is fed into T\ is fed back into input 4 of H. Output 
3 of if is fed into <3>o is fed into T is fed back into input 3 of H. The remaining 
two outputs are reordered." 

To represent the coherently driven dynamics, one then calculates 



N, 



vac 



P (1)0) < [(J 2 ffl (T < $ )) < [(J 3 ffl (Ti < $0) < ff] 4 _J 



(10) 




(11) 
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If one then plugs in a quadrature hybrid model for H (i.e. \x = l/y/%, v = i/y/2) and sets 
both phase shifts to n, the resulting symbolic N = (Sat, Ltv, i/jv) triple is relatively simple 



S 



N 



2y/2i/3 1/3 
1/3 2^/3 

=^a + 2i^ ai + 2i^-a - §e*i 



tj A t At , XO +2 2 i Xl t2 2 , 



/>■ 



a^ao + i- 



'2k 



6 



a Q cti — i- 



/2k Jk , 

-7— ai^o + -s-Oiai + rt(c 
o 3 



where a{o,i} is the annihilation operator for T{ ,i}, analogous labeling applies to Aj and Xj, 
and a{ ,i} are the coherent drive amplitudes driving inputs and 1. 

At this stage, one could produce the the full quantum mechanical equations of motion. 
However, we also invoke a semiclassical approximation that is appropriate for our measure- 
ments in the main Letter. That is, we instead calculate the equations of motion for the 
expectations of the degrees of freedom (e.g. a« = (a*)) and assume that the expectations 



of normal-ordered operators factor (e.g. {a\a,i) 



« \a,.\ 2 : 



. Moreover, as the inputs to N are 



vacuum fields (recall, the coherent drives that excite the network are actually part of N), 
all the noise terms drop out of these expressions and using Eqs. (??-??), we are left with a 
closed system of equations 



d 
dt 
d 
dt 



d 

dt ao 
d 

dt ai 

(A ou tfi) 

\A™t,i) 



(iA + K/3)a - iXoa* al - 2i^^a 1 + -~^(V2ia + a%) 
(iAi + k/ 3)ai - ixi^i - 2i-^r-ao ~ -^-( a o + iVZati] 

o o 

\/2k~ ~-V2 1 

— — a + 2i— a x + 2i— a - -a x 

V2k„ ^2 1 

— — ax + 2i—a - 2%—(Xx + -a . 



(13) 



In the main Letter, the symbolic semiclassical equations of motion analogous to Eqs. (13) 



were produced by QHDL using the slightly more complex schematic in Fig. [7] (which would 
take up pages of complex expressions to reproduce here - symbolic algebra capabilities are 
still a work in progress), which includes transmission line loss and a general phase shift 
parameter. Despite their complexity, when numerical parameters were substituted, the 
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resulting nonlinear, complex equations of motion contained only a small number of terms 
(see below). These equations of motion were typically integrated numerically in minutes on 
a laptop, forming the basis of the simulations presented in the main Letter. 

In the actual model used in the main Letter, the model parameters were // = l/\/2, 
v = phase delays of 2.65 rad, a loss per TKC-pass of 0.4 dB, and the unitless TKC 

parameters k = l/y/3, A = 0.69, and \ — — 4/t 2 /3y / 3 (normalized such that A = 1 and 
a ,i=l correspond to the critical Kerr detuning and drive amplitudes). The equations of 
motion that QHDL produces for these parameters are 

^-a = 0.256600119639834^6^ - 0.2698357229814365 - 0.944502934755685m - 

(JL 

0.117914147124703fii - 0.582406649882899i&i - 0.16747234932605a + 

0.557479734216854m + 0.383245128440802a! - 0.0775918723951391^ 

■f-ai = -0.1179141471247035 - 0.582406649882899ia + 0.2566001 19639834m*^ - 
at 

0.269835722981436^1 - 0.944502934755685^1 - 0.383245128440802a + 
0.077591872395139Ha + 0.16747234932605«i - 0.557479734216854i«i 
^-(Aoutfl) = -0.2861745309459195 + 0.236841667739385ia + 

(JL 

0.109731478271128fii + 0.541990458354401ifii + 0.155850582048067a + 
0.895420212859944ia - 0.356649778754219«i + 0.0722073734777465i«i 
^-(Aout,i) = -0.286174530945919ai + 0.236841667739385mi + 

(JL 

0.1097314782711285 + 0.541990458354401ia + 0.155850582048067«i - 
0.895420212859944iai - 0.356649778754219« - 0.0722073734777465ia - (14) 

Linearized model 

In this section, we primarily describe the linearized model that was used to support the 
hypothesis that transmission line delays are the main cause for the discrepancy between the 
observed and simulated output power oscillations (Fig. 5a & c in the main Letter). We 
thank H. Mabuchi for suggesting the outline of this approach. 

First, though, we give a qualitative argument for the delay-induced discrepancy between 
our system and the model. For low power In drives (and no Ini drive), intra- network signal 
power is too weak for the TKC non-linearity to be significant. According to the QHDL- 
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produced model, as the Ino drive increases past p = 1.2dB, the typical power incident on 
TKC i increases past p = 0.4 (the TKCo incident power is higher still) and the non-linearity 
of both resonators becomes significant (see Fig. |8p,), leading to the sustained oscillations 
encountered in the main Letter. As in the case of bistability, one may roughly understand 
these astable dynamics through a sequence of events. As seen in Fig. (where plot colors 
correspond to the network signals as colored in Fig. [8]d), with a sufficiently strong drive on 
port In , a rise in the signal power incident on TKCo (red) results in an increase in both the 
power exiting Outi and incident on TKCi (green) after a characteristic relaxation time. As 
the TKCi incident power rises past p ~ 0.4, the TKCi reflected signal (purple) begins to 
destructively interfere with the drive signal, causing the TKCo signal to decrease in power 
and the Outo signal (blue) to increase. (Somewhat interestingly, the power of the TKCi 
reflected signal is relatively stable, while its phase - not shown - varies strongly) Eventually, 
as the TKCo incident power drops, the TKCi incident power also drops. The In drive then 
begins to build up the TKC incident power again, and the cycle continues. 

Transmission line delays would allow perturbations to the TKCi incident power to grow 
larger before interferometric feedback is able to counteract them, leading to enhanced insta- 
bility. For example, in the (no-delay) QHDL model, when the Ino drive is set at p = — 1.6dB, 
the steady state TKCi incident power is p = 0.34 (not shown). In the experimental network, 




Time [microseconds] 

FIG. 8: a) Steady state reflected phase from a Kerr resonator driven with the experimentally- 
employed .69A C detuning. Note that the x-axis is in units of power, in contrast to Fig. la in the 
main Letter, b) Depiction of the astable configuration presented in the main Letter, with colored 
arrows indicating various signal segments presented in c. c) Given an Ino drive of p = 2.58dB 
(compare to Fig. 5 in the main Letter) and starting from an un-driven state, a few signal power 
cycles simulated from the QHDL-derived model are plotted. As in figure b, red corresponds to the 
power incident on TKCo, green is both the power exiting Outi and incident on TKCi, purple is 
the signal emitted by TKCi, and blue is the signal exiting Outo- 
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delays add up to 0.48k _1 = 5ns round trip. In the sustained oscillations depicted in Fig. [8p, 
the TKCi incident power increases from p = 0.34 to p — 0.7 in 5ns, well into the non-linear 
regime for the device. Thus, one would expect that transmission line delays would lead to 
stable limit cycles emerging at lower drive powers in general, and that for the experimental 
system at hand, instability with an Ino drive of only p = — 1.6dB would not be unreasonable, 
given typical rates of signal power variation and round trip delays. These expectations are 
given a more quantitative foundation in the remainder of this section. 

The relevant dynamical fixed points of {a,Q,di} for In power drives in the range p = 
{— 2,5}dB were found numerically using Eqs. (14). We then note that for a = Uq + Wq, 



at — Ui + ivi, the equations of motion for {ao, di} from Eqs. (14) may be written as 
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where 1] = 0.256600119639834, A' and B are 4x4 real matrices, 9ft{a:} and are the 

real and imaginary components of a, and we have used ia*af = i{uf + vf){ui + ivi). 
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where we have re-defined the {itj, v j, a{\ now as deviations about the fixed points, and x and 
u are vectors of these deviations. The Ay are 2x2 real matrices, which are dependent on 
the mean In drive through the fixed points. 



Similarly, the definition of the output field fluxes ■jz(A (m t,i) from Eqs. (14) can be written 
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in matrix form as 



y = Cx + Du 



(17) 



where y = [^(v4 onti0 ), J^(A™i,i)] T , and C and D are complex 2x4 matrices. 



We note that the equation of motion Eq. ( 16 ) may be modeled as a linear feedback system 



shown in the dotted box in Fig. |9j This feedback network represents a system whose input 
is the B-transformed input deviations, Bu, and whose output is the deviations of the TKCs' 
internal fields from their fixed points, x. Using the linearized network model represented in 
Fig. [9j we can approximate the consequences of transmission line delays on the overall I/O 
network dynamics about the calculated fixed points (delays will not effect the fixed point 
locations). We make this approximation by inserting delay 'components' e~ ST , where r is 
the time delay and s is the Laplace transform variable, on the feedback lines through which 
the TKCo deviations, Xq, drive the TKCi deviations, x\, and vice versa. This is motivated 
by the intuition that the dominant contributions to these dynamical 'cross terms' have to 
travel 50 cm of SMA cable (two cable interconnections) in order to drive the dynamics 
in the other TKC Note that even within the linearized model this is an approximation. 
For example, the Xq contribution that makes multiple 'passes' through the network before 
driving either x,\ or x are ignored. This approximation is justifiable in that the 'Q' of the 
network is very low - the residual energy left in signals will be low after a few reflections by 
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FIG. 9: Equivalent linearized feedback network suggested by the equations of motion Eqs. ( 16 



Using this linear model, the effects of transmission line delays may be approximated by inserting 
delay 'components' e~ ST in signal lines that represent the driving of the internal state of TKCi by 
the internal state of TKCo and vice versa. 
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the 3 dB hybrid. Using a 5 th -order Pade transfer function approximation of e~ ST and the Ino 
drive-dependent matrices, we can use the Matlab Control Systems Toolbox to calculate 
a minimal state space model for the feedback network depicted in the dotted box in Fig. [9} 
From this model, we can use additional Toolbox functions for I/O pole- zero analysis [36] of 
the entire linearized network depicted in Fig. [9j 

After transforming the linearized dynamics back into dimensionfull parameters, in Fig. 10 
we plot the pole-zero maps for real In drives to Outo signals for two cases: when trans- 
mission line delays are ignored and when they are approximated as described above. The 
other I/O maps produce very similar trends. When no delays are modeled, one sees a 
marginally stable complex conjugate pole pair move towards the imaginary axis as drive 
power increases. At an input drive of p =1.2 dB, this pair crosses the imaginary axis with 
imaginary components ±2ir x 17 MHz, characteristic of a supercritical Hopf bifurcation that 
destabilizes the fixed point to a 17 MHz limit cycle. At higher drives still, the magnitude 
of the imaginary components of this pair keeps increasing, suggesting that the limit cycle 
frequency similarly increases. This interpretation is strongly supported by the simulated 
power spectrum in Fig. 5c in the main Letter: at 1.2 dB drives, 17 MHz power oscillations 
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FIG. 10: I/O pole-zero maps produced by the Matlab Control Systems Toolbox model of the 
linear network represented in Fig. [9] for various Ino drive powers. Poles are marked with x's, zeros 
with o's. Left, when no transmission line delays are included in the linearized model, a complex 
pole pair crosses the imaginary axis at p =1.2 dB drive power and ±2ir x 17 MHz, accurately 
predicting the appearance of the power oscillations observed in model simulation of Fig. 5c in the 
main Letter. Right, when 50 cm transmission line delays are included this pole pair is destabilized, 
crossing the imaginary axis at -1.6 dB drive power and ±2ir x 11 MHz, suggesting the appearance 
of power oscillations much closer to what was observed experimentally in Fig. 5a of the main 
Letter. 
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suddenly appear and increase in frequency with increasing drive power. It is well known that 
the precursors to Hopf bifurcations can be useful for the amplification of AC signals [37J, 
suggesting another potential application for our network. When transmission line delays are 
approximately modeled, the most conspicuous consequence is to further destabilize this pole 
pair. Starting much closer to the imaginary axis, the pair crosses it at -1.6 dB with imag- 
inary components ±2tt x 11 MHz whose magnitudes increase further with increasing drive 
power. This suggests that if transmission line delays were included in the QHDL-produced 
model, 11 MHz power oscillations would first be observed at -1.6 dB in simulation, much 
closer to what was experimentally observed in Fig. 5a of the main Letter. 
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